home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / snpd1292.zip / SUNRISET.C < prev    next >
C/C++ Source or Header  |  1992-12-26  |  22KB  |  508 lines
  1. /*
  2.  
  3. SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
  4.              the length of the day at any date and latitude
  5.  
  6. Written as DAYLEN.C, 1989-08-16
  7.  
  8. Modified to SUNRISET.C, 1992-12-01
  9.  
  10. (c) Paul Schlyter, 1989, 1992
  11.  
  12. Released to the public domain by Paul Schlyter, December 1992
  13.  
  14. */
  15.  
  16.  
  17. #include <stdio.h>
  18. #include <math.h>
  19.  
  20.  
  21. /* A macro to compute the number of days elapsed since 2000 Jan 0.0 */
  22. /* (which is equal to 1999 Dec 31, 0h UT)                           */
  23.  
  24. #define days_since_2000_Jan_0(y,m,d) \
  25.     (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
  26.  
  27. /* Some conversion factors between radians and degrees */
  28.  
  29. #ifndef PI
  30.  #define PI        3.1415926535897932384
  31. #endif
  32.  
  33. #define RADEG     ( 180.0 / PI )
  34. #define DEGRAD    ( PI / 180.0 )
  35.  
  36. /* The trigonometric functions in degrees */
  37.  
  38. #define sind(x)  sin((x)*DEGRAD)
  39. #define cosd(x)  cos((x)*DEGRAD)
  40. #define tand(x)  tan((x)*DEGRAD)
  41.  
  42. #define atand(x)    (RADEG*atan(x))
  43. #define asind(x)    (RADEG*asin(x))
  44. #define acosd(x)    (RADEG*acos(x))
  45. #define atan2d(y,x) (RADEG*atan2(y,x))
  46.  
  47.  
  48. /* Following are some macros around the "workhorse" function __daylen__ */
  49. /* They mainly fill in the desired values for the reference altitude    */
  50. /* below the horizon, and also selects whether this altitude should     */
  51. /* refer to the Sun's center or its upper limb.                         */
  52.  
  53.  
  54. /* This macro computes the length of the day, from sunrise to sunset. */
  55. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  56. /* 35 arc minutes below the horizon (this accounts for the refraction */
  57. /* of the Earth's atmosphere).                                        */
  58. #define day_length(year,month,day,lon,lat)  \
  59.         __daylen__( year, month, day, lon, lat, -35.0/60.0, 1 )
  60.  
  61. /* This macro computes the length of the day, including civil twilight. */
  62. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  63. /* the horizon.                                                         */
  64. #define day_civil_twilight_length(year,month,day,lon,lat)  \
  65.         __daylen__( year, month, day, lon, lat, -6.0, 0 )
  66.  
  67. /* This macro computes the length of the day, incl. nautical twilight.  */
  68. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  69. /* below the horizon.                                                   */
  70. #define day_nautical_twilight_length(year,month,day,lon,lat)  \
  71.         __daylen__( year, month, day, lon, lat, -12.0, 0 )
  72.  
  73. /* This macro computes the length of the day, incl. astronomical twilight. */
  74. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  75. /* below the horizon.                                                      */
  76. #define day_astronomical_twilight_length(year,month,day,lon,lat)  \
  77.         __daylen__( year, month, day, lon, lat, -18.0, 0 )
  78.  
  79.  
  80. /* This macro computes times for sunrise/sunset.                      */
  81. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  82. /* 35 arc minutes below the horizon (this accounts for the refraction */
  83. /* of the Earth's atmosphere).                                        */
  84. #define sun_rise_set(year,month,day,lon,lat,rise,set)  \
  85.         __sunriset__( year, month, day, lon, lat, -35.0/60.0, 1, rise, set )
  86.  
  87. /* This macro computes the start and end times of civil twilight.       */
  88. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  89. /* the horizon.                                                         */
  90. #define civil_twilight(year,month,day,lon,lat,start,end)  \
  91.         __sunriset__( year, month, day, lon, lat, -6.0, 0, start, end )
  92.  
  93. /* This macro computes the start and end times of nautical twilight.    */
  94. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  95. /* below the horizon.                                                   */
  96. #define nautical_twilight(year,month,day,lon,lat,start,end)  \
  97.         __sunriset__( year, month, day, lon, lat, -12.0, 0, start, end )
  98.  
  99. /* This macro computes the start and end times of astronomical twilight.   */
  100. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  101. /* below the horizon.                                                      */
  102. #define astronomical_twilight(year,month,day,lon,lat,start,end)  \
  103.         __sunriset__( year, month, day, lon, lat, -18.0, 0, start, end )
  104.  
  105.  
  106. /* Function prototypes */
  107.  
  108. double __daylen__( int year, int month, int day, double lon, double lat,
  109.                    double altit, int upper_limb );
  110.  
  111. int __sunriset__( int year, int month, int day, double lon, double lat,
  112.                   double altit, int upper_limb, double *rise, double *set );
  113.  
  114. void sunpos( double d, double *lon, double *r );
  115.  
  116. void sun_RA_dec( double d, double *RA, double *dec, double *r );
  117.  
  118. double revolution( double x );
  119.  
  120. double rev180( double x );
  121.  
  122. double GMST0( double d );
  123.  
  124.  
  125.  
  126. /* A small test program */
  127.  
  128. void main(void)
  129. {
  130.       int year,month,day;
  131.       double lon, lat;
  132.       double daylen, civlen, nautlen, astrlen;
  133.       double rise, set, civ_start, civ_end, naut_start, naut_end,
  134.              astr_start, astr_end;
  135.       int    rs, civ, naut, astr;
  136.  
  137.       printf( "Longitude (+ is east) and latitude (+ is north) : " );
  138.       scanf( "%lf %lf", &lon, &lat );
  139.  
  140.       for(;;)
  141.       {
  142.             printf( "Input date ( yyyy mm dd ) (ctrl-C exits): " );
  143.             scanf( "%d %d %d", &year, &month, &day );
  144.  
  145.             daylen  = day_length(year,month,day,lon,lat);
  146.             civlen  = day_civil_twilight_length(year,month,day,lon,lat);
  147.             nautlen = day_nautical_twilight_length(year,month,day,lon,lat);
  148.             astrlen = day_astronomical_twilight_length(year,month,day,
  149.                   lon,lat);
  150.  
  151.             printf( "Day length:                 %5.2f hours\n", daylen );
  152.             printf( "With civil twilight         %5.2f hours\n", civlen );
  153.             printf( "With nautical twilight      %5.2f hours\n", nautlen );
  154.             printf( "With astronomical twilight  %5.2f hours\n", astrlen );
  155.             printf( "Length of twilight: civil   %5.2f hours\n",
  156.                   (civlen-daylen)/2.0);
  157.             printf( "                  nautical  %5.2f hours\n",
  158.                   (nautlen-daylen)/2.0);
  159.             printf( "              astronomical  %5.2f hours\n",
  160.                   (astrlen-daylen)/2.0);
  161.  
  162.             rs   = sun_rise_set         ( year, month, day, lon, lat,
  163.                                           &rise, &set );
  164.             civ  = civil_twilight       ( year, month, day, lon, lat,
  165.                                           &civ_start, &civ_end );
  166.             naut = nautical_twilight    ( year, month, day, lon, lat,
  167.                                           &naut_start, &naut_end );
  168.             astr = astronomical_twilight( year, month, day, lon, lat,
  169.                                           &astr_start, &astr_end );
  170.  
  171.             printf( "Sun at south %5.2fh UT\n", (rise+set)/2.0 );
  172.  
  173.             switch( rs )
  174.             {
  175.                 case 0:
  176.                     printf( "Sun rises %5.2fh UT, sets %5.2fh UT\n",
  177.                              rise, set );
  178.                     break;
  179.                 case +1:
  180.                     printf( "Sun above horizon\n" );
  181.                     break;
  182.                 case -1:
  183.                     printf( "Sun below horizon\n" );
  184.                     break;
  185.             }
  186.  
  187.             switch( civ )
  188.             {
  189.                 case 0:
  190.                     printf( "Civil twilight starts %5.2fh, "
  191.                             "ends %5.2fh UT\n", civ_start, civ_end );
  192.                     break;
  193.                 case +1:
  194.                     printf( "Never darker than civil twilight\n" );
  195.                     break;
  196.                 case -1:
  197.                     printf( "Never as bright as civil twilight\n" );
  198.                     break;
  199.             }
  200.  
  201.             switch( naut )
  202.             {
  203.                 case 0:
  204.                     printf( "Nautical twilight starts %5.2fh, "
  205.                             "ends %5.2fh UT\n", naut_start, naut_end );
  206.                     break;
  207.                 case +1:
  208.                     printf( "Never darker than nautical twilight\n" );
  209.                     break;
  210.                 case -1:
  211.                     printf( "Never as bright as nautical twilight\n" );
  212.                     break;
  213.             }
  214.  
  215.             switch( astr )
  216.             {
  217.                 case 0:
  218.                     printf( "Astronomical twilight starts %5.2fh, "
  219.                             "ends %5.2fh UT\n", astr_start, astr_end );
  220.                     break;
  221.                 case +1:
  222.                     printf( "Never darker than astronomical twilight\n" );
  223.                     break;
  224.                 case -1:
  225.                     printf( "Never as bright as astronomical twilight\n" );
  226.                     break;
  227.             }
  228.  
  229.       }
  230. }
  231.  
  232.  
  233. /* The "workhorse" function for sun rise/set times */
  234.  
  235. int __sunriset__( int year, int month, int day, double lon, double lat,
  236.                   double altit, int upper_limb, double *trise, double *tset )
  237. /***************************************************************************/
  238. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  239. /*       Eastern longitude positive, Western longitude negative       */
  240. /*       Northern latitude positive, Southern latitude negative       */
  241. /*       The longitude value IS critical in this function!            */
  242. /*       altit = the altitude which the Sun should cross              */
  243. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  244. /*               for civil, -12 degrees for nautical and -18          */
  245. /*               degrees for astronomical twilight.                   */
  246. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  247. /*               Set to non-zero (e.g. 1) when computing rise/set     */
  248. /*               times, and to zero when computing start/end of       */
  249. /*               twilight.                                            */
  250. /*        *rise = where to store the rise time                        */
  251. /*        *set  = where to store the set  time                        */
  252. /*                Both times are relative to the specified altitude,  */
  253. /*                and thus this function can be used to comupte       */
  254. /*                various twilight times, as well as rise/set times   */
  255. /* Return value:  0 = sun rises/sets this day, times stored at        */
  256. /*                    *trise and *tset.                               */
  257. /*               +1 = sun above the specified "horizon" 24 hours.     */
  258. /*                    *trise set to time when the sun is at south,    */
  259. /*                    minus 12 hours while *tset is set to the south  */
  260. /*                    time plus 12 hours. "Day" length = 24 hours     */
  261. /*               -1 = sun is below the specified "horizon" 24 hours   */
  262. /*                    "Day" length = 0 hours, *trise and *tset are    */
  263. /*                    both set to the time when the sun is at south.  */
  264. /*                                                                    */
  265. /**********************************************************************/
  266. {
  267.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  268.       sr,         /* Solar distance, astronomical units */
  269.       sRA,        /* Sun's Right Ascension */
  270.       sdec,       /* Sun's declination */
  271.       sradius,    /* Sun's apparent radius */
  272.       t,          /* Diurnal arc */
  273.       tsouth,     /* Time when Sun is at south */
  274.       sidtime;    /* Local sidereal time */
  275.  
  276.       int rc = 0; /* Return cde from function - usually 0 */
  277.  
  278.       /* Compute d of 12h local mean solar time */
  279.       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  280.  
  281.       /* Compute local sideral time of this moment */
  282.       sidtime = revolution( GMST0(d) + 180.0 + lon );
  283.  
  284.       /* Compute Sun's RA + Decl at this moment */
  285.       sun_RA_dec( d, &sRA, &sdec, &sr );
  286.  
  287.       /* Compute time when Sun is at south - in hours UT */
  288.       tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
  289.  
  290.       /* Compute the Sun's apparent radius, degrees */
  291.       sradius = 0.2666 / sr;
  292.  
  293.       /* Do correction to upper limb, if necessary */
  294.       if ( upper_limb )
  295.             altit -= sradius;
  296.  
  297.       /* Compute the diurnal arc that the Sun traverses to reach */
  298.       /* the specified altitide altit: */
  299.       {
  300.             double cost;
  301.             cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
  302.                   ( cosd(lat) * cosd(sdec) );
  303.             if ( cost >= 1.0 )
  304.                   rc = -1, t = 0.0;       /* Sun always below altit */
  305.             else if ( cost <= -1.0 )
  306.                   rc = +1, t = 12.0;      /* Sun always above altit */
  307.             else
  308.                   t = acosd(cost)/15.0;   /* The diurnal arc, hours */
  309.       }
  310.  
  311.       /* Store rise and set times - in hours UT */
  312.       *trise = tsouth - t;
  313.       *tset  = tsouth + t;
  314.  
  315.       return rc;
  316. }  /* __sunriset__ */
  317.  
  318.  
  319.  
  320. /* The "workhorse" function */
  321.  
  322.  
  323. double __daylen__( int year, int month, int day, double lon, double lat,
  324.                    double altit, int upper_limb )
  325. /**********************************************************************/
  326. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  327. /*       Eastern longitude positive, Western longitude negative       */
  328. /*       Northern latitude positive, Southern latitude negative       */
  329. /*       The longitude value is not critical. Set it to the correct   */
  330. /*       longitude if you're picky, otherwise set to to, say, 0.0     */
  331. /*       The latitude however IS critical - be sure to get it correct */
  332. /*       altit = the altitude which the Sun should cross              */
  333. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  334. /*               for civil, -12 degrees for nautical and -18          */
  335. /*               degrees for astronomical twilight.                   */
  336. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  337. /*               Set to non-zero (e.g. 1) when computing day length   */
  338. /*               and to zero when computing day+twilight length.      */
  339. /**********************************************************************/
  340. {
  341.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  342.       obl_ecl,    /* Obliquity (inclination) of Earth's axis */
  343.       sr,         /* Solar distance, astronomical units */
  344.       slon,       /* True solar longitude */
  345.       sin_sdecl,  /* Sine of Sun's declination */
  346.       cos_sdecl,  /* Cosine of Sun's declination */
  347.       sradius,    /* Sun's apparent radius */
  348.       t;          /* Diurnal arc */
  349.  
  350.       /* Compute d of 12h local mean solar time */
  351.       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  352.  
  353.       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
  354.       obl_ecl = 23.4393 - 3.563E-7 * d;
  355.  
  356.       /* Compute Sun's position */
  357.       sunpos( d, &slon, &sr );
  358.  
  359.       /* Compute sine and cosine of Sun's declination */
  360.       sin_sdecl = sind(obl_ecl) * sind(slon);
  361.       cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
  362.  
  363.       /* Compute the Sun's apparent radius, degrees */
  364.       sradius = 0.2666 / sr;
  365.  
  366.       /* Do correction to upper limb, if necessary */
  367.       if ( upper_limb )
  368.             altit -= sradius;
  369.  
  370.       /* Compute the diurnal arc that the Sun traverses to reach */
  371.       /* the specified altitide altit: */
  372.       {
  373.             double cost;
  374.             cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
  375.                   ( cosd(lat) * cos_sdecl );
  376.             if ( cost >= 1.0 )
  377.                   t = 0.0;                      /* Sun always below altit */
  378.             else if ( cost <= -1.0 )
  379.                   t = 24.0;                     /* Sun always above altit */
  380.             else  t = (2.0/15.0) * acosd(cost); /* The diurnal arc, hours */
  381.       }
  382.       return t;
  383. }  /* __daylen__ */
  384.  
  385.  
  386. /* This function computes the Sun's position at any instant */
  387.  
  388. void sunpos( double d, double *lon, double *r )
  389. /******************************************************/
  390. /* Computes the Sun's ecliptic longitude and distance */
  391. /* at an instant given in d, number of days since     */
  392. /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
  393. /* computed, since it's always very near 0.           */
  394. /******************************************************/
  395. {
  396.       double M,         /* Mean anomaly of the Sun */
  397.              w,         /* Mean longitude of perihelion */
  398.                         /* Note: Sun's mean longitude = M + w */
  399.              e,         /* Eccentricity of Earth's orbit */
  400.              E,         /* Eccentric anomaly */
  401.              x, y,      /* x, y coordinates in orbit */
  402.              v;         /* True anomaly */
  403.  
  404.       /* Compute mean elements */
  405.       M = revolution( 356.0470 + 0.9856002585 * d );
  406.       w = 282.9404 + 4.70935E-5 * d;
  407.       e = 0.016709 - 1.151E-9 * d;
  408.  
  409.       /* Compute true longitude and radius vector */
  410.       E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
  411.             x = cosd(E) - e;
  412.       y = sqrt( 1.0 - e*e ) * sind(E);
  413.       *r = sqrt( x*x + y*y );              /* Solar distance */
  414.       v = atan2d( y, x );                  /* True anomaly */
  415.       *lon = v + w;                        /* True solar longitude */
  416.       if ( *lon >= 360.0 )
  417.             *lon -= 360.0;                   /* Make it 0..360 degrees */
  418. }
  419.  
  420. void sun_RA_dec( double d, double *RA, double *dec, double *r )
  421. {
  422.       double lon, obl_ecl, x, y, z;
  423.  
  424.       /* Compute Sun's ecliptical coordinates */
  425.       sunpos( d, &lon, r );
  426.  
  427.       /* Compute ecliptic rectangular coordinates (z=0) */
  428.       x = *r * cosd(lon);
  429.       y = *r * sind(lon);
  430.  
  431.       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
  432.       obl_ecl = 23.4393 - 3.563E-7 * d;
  433.  
  434.       /* Convert to equatorial rectangular coordinates - x is uchanged */
  435.       z = y * sind(obl_ecl);
  436.       y = y * cosd(obl_ecl);
  437.  
  438.       /* Convert to spherical coordinates */
  439.       *RA = atan2d( y, x );
  440.       *dec = atan2d( z, sqrt(x*x + y*y) );
  441.  
  442. }  /* sun_RA_dec */
  443.  
  444.  
  445. /******************************************************************/
  446. /* This function reduces any angle to within the first revolution */
  447. /* by subtracting or adding even multiples of 360.0 until the     */
  448. /* result is >= 0.0 and < 360.0                                   */
  449. /******************************************************************/
  450.  
  451. #define INV360    ( 1.0 / 360.0 )
  452.  
  453. double revolution( double x )
  454. /*****************************************/
  455. /* Reduce angle to within 0..360 degrees */
  456. /*****************************************/
  457. {
  458.       return( x - 360.0 * floor( x * INV360 ) );
  459. }  /* revolution */
  460.  
  461. double rev180( double x )
  462. /*********************************************/
  463. /* Reduce angle to within +180..+180 degrees */
  464. /*********************************************/
  465. {
  466.       return( x - 360.0 * floor( x * INV360 + 0.5 ) );
  467. }  /* revolution */
  468.  
  469.  
  470. /*******************************************************************/
  471. /* This function computes GMST0, the Greenwhich Mean Sidereal Time */
  472. /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  */
  473. /* 0h UT).  GMST is then the sidereal time at Greenwich at any     */
  474. /* time of the day.  I've generelized GMST0 as well, and define it */
  475. /* as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at */
  476. /* other times than 0h UT as well.  While this sounds somewhat     */
  477. /* contradictory, it is very practical:  instead of computing      */
  478. /* GMST like:                                                      */
  479. /*                                                                 */
  480. /*  GMST = (GMST0) + UT * (366.2422/365.2422)                      */
  481. /*                                                                 */
  482. /* where (GMST0) is the GMST last time UT was 0 hours, one simply  */
  483. /* computes:                                                       */
  484. /*                                                                 */
  485. /*  GMST = GMST0 + UT                                              */
  486. /*                                                                 */
  487. /* where GMST0 is the GMST "at 0h UT" but at the current moment!   */
  488. /* Defined in this way, GMST0 will increase with about 4 min a     */
  489. /* day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)   */
  490. /* is equal to the Sun's mean longitude plus/minus 180 degrees!    */
  491. /* (if we neglect aberration, which amounts to 20 seconds of arc   */
  492. /* or 1.33 seconds of time)                                        */
  493. /*                                                                 */
  494. /*******************************************************************/
  495.  
  496. double GMST0( double d )
  497. {
  498.       double sidtim0;
  499.       /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr  */
  500.       /* L = M + w, as defined in sunpos().  Since I'm too lazy to */
  501.       /* add these numbers, I'll let the C compiler do it for me.  */
  502.       /* Any decent C compiler will add the constants at compile   */
  503.       /* time, imposing no runtime or code overhead.               */
  504.       sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) +
  505.                           ( 0.9856002585 + 4.70935E-5 ) * d );
  506.       return sidtim0;
  507. }  /* GMST0 */
  508.